package weka.em;

import static java.lang.Math.exp;
import static java.lang.Math.log;

public class LDAUtils
{
	public static double digamma(double x)
	{
		double p;
		assert x > 0;
		x = x + 6;
		p = 1 / (x * x);
		p = (((0.004166666666667 * p - 0.003968253986254) * p + 0.008333333333333)
				* p - 0.083333333333333)
				* p;
		p = p + log(x) - 0.5 / x - 1 / (x - 1) - 1 / (x - 2) - 1 / (x - 3) - 1
				/ (x - 4) - 1 / (x - 5) - 1 / (x - 6);
		return p;
	}

	public static double trigamma(double x)
	{
		double result;
		double neginf = -1.0 / 0.0, small = 1e-4, large = 8, c = 1.6449340668482264365,
		/*
		 * pi^ 2 / 6 = Zeta ( 2 )
		 */
		c1 = -2.404113806319188570799476, /* -2 Zeta(3) */
		b2 = 1. / 6, b4 = -1. / 30, b6 = 1. / 42, b8 = -1. / 30, b10 = 5. / 66;

		/* Negative values */
		/*
		 * Use the derivative of the digamma reflection formula: -trigamma(-x) =
		 * trigamma(x+1) - (pi*csc(pi*x))^2
		 */
		if (x < 0)
		{
			result = Math.PI / Math.sin(-Math.PI * x);
			return -trigamma(1 - x) + result * result;
		}
		/* Use Taylor series if argument <= small */
		if (x <= small)
		{
			return 1 / (x * x) + c + c1 * x;
		}
		result = 0;
		/* Reduce to trigamma(x+n) where ( X + N ) >= B */
		while (x < large)
		{
			result += 1 / (x * x);
			x++;
		}
		/* Apply asymptotic formula when X >= B */
		/* This expansion can be computed in Maple via asympt(Psi(1,x),x) */
		if (x >= large)
		{
			double r = 1 / (x * x);
			result += 0.5 * r
					+ (1 + r * (b2 + r * (b4 + r * (b6 + r * (b8 + r * b10)))))
					/ x;
		}
		return result;
	}

	static boolean GAMMAACCURATE = true;

	static final double HALFLN2PI = 9.189385332046727E-001;

	static final double HALFLN2 = 3.465735902799726e-001;

	static final double INV810 = 1.234567901234568E-003;

	public static double lgamma(double x)
	{
		double lnx, einvx;
		double prec;
		assert x > 0;
		lnx = log(x);
		einvx = exp(1. / x);

		if (GAMMAACCURATE)
		{
			prec = x * x * x;
			prec *= prec;
			prec = INV810 / prec;
			/*
			 * y = x * ( log(x) - 1 + .5 * log(x * sinh(1/x) + prec) ) - .5 *
			 * log(x) + .5 * log(2 * pi)
			 */
			return x
					* (lnx - 1. + .5 * log(x * (einvx - 1. / einvx) / 2. + prec))
					- .5 * lnx + HALFLN2PI;
		}
		else
		{
			/*
			 * y = x * ( 1.5 * log(x) - 1 + .5 * log(exp(1/x) - 1/exp(1/x)) - .5
			 * * log(2) ) - .5 * log(x) + .5 * log(2 * pi);
			 */
			return x
					* (1.5 * lnx - 1. + .5 * log(einvx - 1. / einvx) - HALFLN2)
					- .5 * lnx + HALFLN2PI;
		}
	}

	public static double log_gamma(double x)
	{
		double z;
		assert x > 0;
		z = 1. / (x * x);

		x = x + 6;
		z = (((-0.000595238095238 * z + 0.000793650793651) * z - 0.002777777777778)
				* z + 0.083333333333333)
				/ x;
		z = (x - 0.5) * log(x) - x + 0.918938533204673 + z - log(x - 1)
				- log(x - 2) - log(x - 3) - log(x - 4) - log(x - 5)
				- log(x - 6);
		return z;
	}
}
